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Most  variance  reduction  techniques  encountered  in  simulations  are 
designed  to  operate  with  stationary  models.  However,  many  simulations 
are  nonstationary  In  character,  population  growth  mooels  being  an  example. 
One  way  to  facilitate  statistical  inference  in  a nonstationary  simulation 
is  to  interchange  the  order  of  replication  collection  and  time  evolution. 
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That  1s,  at  each  time  point  several  replications  are  performed  to  enable  a 
user  to  estiniate  parameters  at  that  point  in  time.  This  paper  describes  three 
variance  reduction  techniques  that  use  this  interchange  between  collection 
and  evolution  to  induce  negative  correlation  between  replications,  thereby 
producing  estimates  with  smaller  variances.  Model  1 describes  a procedure 
that  occasionally  relies  on  the  solution  of  a linear  program  to  develop  an 
optimal  sampling  plan.  Model  2 offers  an  alternative  that  applies  when  the 
populations  in  strata  are  large.  Model  3 applies  when  survival  probabilities 
are  functions  of  an  exoganous  random  variable  such  as  rainfall.  A female 
elephant  population  simulation  Illustrates  the  success  one  can  expect  with 
morel  1 . 
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Abstract 


Most  variance  reduction  techniques  encountered  in  simulations 
are  designed  to  operate  with  stationary  models.  However,  many  simula- 
tions are  nonstationary  in  character,  population  growth  models  being 
an  e.''.3snple.  One  way  to  facilitate  statistical  inference  in  a nonsta- 
tionary simulation  is  to  interchante  the  order  of  replication  collection 
and  time  evolution.  That  is,  at  each  time  point  several  replications 
are  performed  to  enable  a user  to  estimate  parameters  at  that  point  in 
time.  This  paper  describes  three  variance  reduction  techniques  that  use 
this  interchange  between  collection  and  evolution  to  induce  negative  cor- 
relation between  replications,  thereby  producing  estimates  with  smaller 
variances.  Model  1 describes  a procedure  that  occasionally  relies  on  the 
solution  of  a linear  program  to  develop  an  optimal  sampling  plan.  Model 
2 offers  an  alternative  that  applies  when  the  populations  in  strata  are 
large.  Model  3 applies  when  survival  probabilities  are  functions  of  an 
exogenous  random  variable  such  as  rainfall.  A fe.nale  elephant  population 
simulation  illustrates  the  success  one  can  expect  with  model  1. 


1 . Introduction 


A^iong  the  tools  available  to  the  simulation  user,  variance  reductio/7 
techniques  have  from  their  inception  held  considerable  attraction.  Broao  y 
speaking  a variance  reduction  technique  is  a sampling  plan  that  enables  a user 
to  achieve  a specified  level  of  statistical  accuracy  for  less  cost  than  pure 
random  sampling  permits  or,  conversely,  it  allows  a user  to  realize  greater 
accuracy  for  specified  cost.  Hammersley  and  Handscomb  [2]  describe  variance 
reduction  techniques  applicable  in  Honte  Carlo  simulation.  Fishman  [1], 
Kleijnen  [3]  and  Mihram  [4]  describe  the  extension  of  these  methods  to  dis- 
crete event  digital  simulation. 

In  reviewing  the  techniques  and  the  problems  for  which  they  are  appro- 
priate one  notes  an  emphasis  on  their  use  in  simulations  of  stationary 
models.  A stationary  model  describes  an  environment  in  which  means  and 
covariance  structures  remain  constant  in  time.  However,  the  simulation 
meth»->d  also  applies  to  nonstationary  systems  of  which  population  growth 
models  form  a substantive  share.  In  a population  model  each  year  or  period 
begins  with  a frequency  distribution  of  population  by  stratum  based  on  age 
sex,  etc.  In  advancing  time  by  one  period  the  simulation  exposes  each 
member  of  each  stratum  to  the  risk  of  death  so  that  at  the  ind  of  the  period 
a new  distribution  by  strata  appears.  If  the  expected  frequency  in  each 
stratum  is  constant  through  time,  then  stationarity  prevails.  However,  this 
stationarity  requires  a fortuitous  balance  oetween  birth  rates,  age  specific 

^This  work  was  motivated  ^y  discussions  on  the  construction  of  population 
simulations  for  long-lived  species  with  Professor  Daniel  Botkin  of  the 
Ecosystems  Center  of  the  Marine  Biological  Laboratory  at  Woods  Hole,  Massa- 
chusetts and  Professors  Richard  Miller  and  Matthew  Sobel  of  Yale  University. 

I am  grateful  to  thefr;  for  introducing  me  to  the  problem  and  for  their 
valuable  ir,  .ights.  I am  also  grateful  to  William  Kwapil  and  Mark  Miller 
for  providing  computaf  ’nal  assistance  in  preparing  Tables  1 and  2. 


survival  rates  and  other  rates  More  commonly,  the  shape  of  the  strata 
frequency  distribution,  the  mean  number  of  members  in  the  population  or 
both  change  time  evolves  so  that  nonstationarity  prevai’.s. 

Now  the  difficulty  with  nonstationarity  in  practice  is  that  one  can 
not  apply  the  statistical  metiods  developed  for  analy<:1ng  simulation  output 
for  stationary  models.  These  latter  techniques  take  two  principal  forms. 

The  first  analyzes  data  within  a single  run  to  estimate  parameters  of  interest 
that  are  time  independent  and  the  second  uses  averages  formed  on  individual 
replications  to  estimate  the  parameters  [1].  Because  nons\".ationarity  implies 
time  dependent  pare*»'eters  one  has  to  look  elsewhere  for  analysis  techniques. 
One  approach  is  to  use  the  data  on  sample  frequency  distributions  collected 
on  several  replications  at  the  sai.v*'  time  point  to  estimate  the  true  fre- 
quency distribution  at  that  time.  This  technique  is  legitimate,  provided 
that  all  replications  begin  with  the  same  initia’  strata  frequer.cy  distri- 
bution and  population. 

Analysis  across  replications  can  proceed  in  at  least  two  ways.  One 
way  has  the  simulation  user  run  replications  serially,  store  the  sample 
frequency  data  at  each  time  point  in  each  run  and  then  combine  correspond- 
ing data,  after  all  runs  are  completed,  to  estimate  the  mean  rtrata 
frequ  cy  distribution  at  each  point.  For  a model  with  many  strata  run  for 
a long  . 'riod  of  time  the  raw  data  may  require  substantial  storage.  Moreover, 
the  simulation  user  must  wait  f ,r  co(.iplet1on  of  all  runs  before  learning 
about  interim  estimates. 

In  an  alternative  method  of  analysis  that  we  plan  to  use  here,  the 
simulation  user  perform*  replications  at  a given  time  point  before  proceed- 
ing to  the  next  time.  In  effect,  he  reverses  the  order  of  time  advance  and 
replication.  Doing  so  nables  him  to  estimate  the  mean  strata  frequency 
distribution  at  a given  time  point  before  moving  forward,  if  an  online 
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graphics  device  is  available,  looking  at  successive  iinages  of  this  estimated 
distribution  may  provide  sufficient  Information  to  terminate  the  simulation 
prematurely.  Or  the  information  gained  may  encourage  the  user  to  change 
certain  p.irameters  of  the  simulation  to  see  how  the  population  responds. 

Although  the  display  and  interpretive  feature  alone  make  this  second 
alternative  attractive,  the  ability  to  apply  a variance  reduction  te  nnique 
to  corresponding  strata  in  the  different  replications  at  a given  time  point 
offer  an  additional  attraction.  Especially  since  the  accuracy  achieved 
for  a given  cost  can  be  considerably  greater  than  independent  replications 
allow. 

In  this  paper  we  describe  three  methods  of  variance  reduction  applicable 
to  population  growth  models.  Kodel  1 in  Section  2 outlines  a method  appro- 
priate for  a simulation  in  whic!‘  the  probability  of  survival  in  each  strata 
in  each  year  randomly  assumes  ont  of  two  possible  values.  As  n the 
number  of  replications  increases,  the  achievable  variance  reduction  often 
increases  dramatically.  Occasionally  the  technique  calls  for  the  solution 
of  a moderate  size  linear  program,  perhaps  an  impediment  today  but  not 
necessarily  so  tomorrow. 

Section  3 describes  a second  variance  reduction  technique  that  applies 
for  the  survival  probability  mechanism  of  model  1 when  the  number  in  each 
stratum  is  large.  This  technique  appears  most  relevant  when  the  two  possible 
survival  probabilities  per  strata  are  close  in  value.  Section  4 offers  a 
third  variance  reduction  technique  that  applies  when  the  survival 
probability  is  a continuous  function  of  an  exogenous  variable  such  as  rain- 
fall which  itself  follows  a normal  probability  law.  To  give  the  reader  an 
appreciation  of  potential  perfom3dnce  Section  5 describes  the  application  of 
model  1 in  Section  3 to  a simulation  of  female  elephant  population  growth. 

The  results  are  encouraging. 
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n 


2.  Model  1 


Let  W denote  a random  /ariable  whose  distribution  is  rectangular  on 
(0,1).  Then  W defines  a uniform  deviate  and  U(0,1),  its  distribution.  Let 
B(N,o)  denote  the  binomial  distribution  with  integral  parameter  N and  probability 
0.  Define  X and  Y as  random  variables  from  8(N,o)  and  B(N,B),  respectively, 
and  let  I denote  the  binary  random  variable 


(1) 


1 


I = 


p < W < 1 


To  illustrate  the  use  of  these  concepts,  assume  that  in  a given  year  a 
member  of  a particular  age  stratum  of  a population  has  a survival  probability  a 
with  probability  p or  g with  probability  1-p.  For  example,  in  an  animal 
population  the  different  survival  rates  arise  as  a consequence  of  favorable  (o) 
or  unfavorable  (6)  weather  conditions.  Then  one  expects  a>$  . If  the  stratum 
has  N members  then  the  number  who  survive  to  enter  the  next  age  stratum  the 
next  year  is 

(2)  R = IX  + (l-I)Y  = I(X-Y)  + Y . 


Define  the  operator  a to  mean  deviation  from  expectation.  Then 


(3) 


aR  = N(a-g)al  + p (aX  -aY)  + al(aX-aY)  + aY 


which,  provided  W is  independent  of  X and  Y,  has  mean  zero  and  variance 


(4) 


V = var  (R)  = p(l-p)  (a-g)"^  + N[o(l-o)  + (l-plgCl-e)]. 


Applying  the  procedure  to  each  stratum  with  its  corresponding  p,  N,  a and  e, 
one  can  construct  a sampl e age  distribution  for  the  year  following  the  one  m 
which  sampling  occurred.  Repeated  application  through  successive  years  allows 
one  to  see  how  the  sample  age  distributions  change  over  time. 

Since  our  principal  interest  concerns  nonstationary  features  of  the 
distribution,  it  is  important  that  the  variance  of  the  ordinates  in  each 
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sample  age  distribution  in  each  year  regain  within  acceptable  limits  for  meaninii'- 
ful  inference.  If  in  each  year  one  performs  r.  independent  replications  for 
each  age  stratum,  then  the  sample  mean  ordinate  for  a giver,  stratum  has  variance 
V/n,  V being  determined  by  the  stratum.  The  objective  of  the  present  research 
is  to  devise  a sampling  plan  among  the  n replications  tnat  not  only  produces  a 
smaller  variance  but  actually  produces  the  smallest  possible  for  n replications 
Although  an  attempt  at  variance  reduction  is  always  desirable,  it  assumes  par- 
ticular importance  in  the  present  case.  Notice  that  the  leading  term  in  (4) 
is  proportional  to  N , not  N as  would  occur  in  the  binomial  case  This 

source  of  variation  provides  additional  impetus  in  the  search  for  .s  sampling  plan 
that  will  lead  to  acceptable  estimates  of  the  ordinates  of  an  distribution. 

Let  a subscript  on  I,  N,  X,  Y,  V and  W denote  replication  number. 

Define  W^.  as 

(5)  W.  = 0.+  U (mod  1) 

\ 

0 £ 0^  < 1 i = 1 ,. . . ,n 

so  that  W^,...,  are  uniform,  but  not  independent,  deviates.  Then  for  a 
particular  age  stratum  replication  i yields 


(6) 

(7) 


aR 


Rj  = !,  (X,  - V,)  * Y, 

, = Nj  {o-6)aI,  ♦ p{4X--sYj)  + 4lj  (iXj-4V.)  ♦ 4Yj 


Here  denotes  the  number  of  survivors  on  replication  i and 

— -1 

R^  = n • £ R. 

" i=l  ^ 

denotes  the  sample  mean  number  of  survivors  over  the  n replications. 
If  these  replications  are  independent  then  R^  has  variance  V/n. 


-6- 


Consider  the  case  in  which  X,,...,  X are  i.i.d.,  Y, ,...,  Y are  i.i.d. 

* n in 

and  {X.}  and  {Y.}  are  independent.  Then  R has  variance 
■ I n 

var^^^  (R„)  = {vJ^^+2U^^h/n 


(8) 

(9) 

„v(') 

00) 

nuj” 

' ' = p(l-p)  (a-S)"  Z N-  + [pa(l-a)  + (1-p)  e(l-e)]  Z N- 
" i=l  ’ i=l  ' 

y n-1  n 

nu;"  = (a-B)^  £ Z N.  N-  cov  (I.,T.) 

" i=l  j=i+l  ^ ^ ^ ^ 

where  the  superscript^^ ^denotes  the  model  number.  Let 

i max  (0,4.)  . 

Since  a little  thought  shows  that  for  {e^.  i=2,...,  n) 

OU  E(I-  Ij.)  = (p-ej''^  + (p-l-e.+0j)^  + (p+e,-0j-(p-9j)'^)'^  j>i, 
one  can  show  that  the  optimal  0^...,  for  minimizing  U^- ^emerge 


from  the  solution  of  the  linear  program: 

n-1  n 


(12) 


= min  E E N.  N-  (D-  + E-. 
9 i=l  j=i+l  ^ J J 


1 (Op...,  0^) 


subject  to  the  n(n+l)  constraints 


9]  1 0 


(13) 


°j  - “j-1 


“j  i p-“j  j 
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* j = i+1 ,. . . , n;  i=l , . . . , n-1 

F,  . . P«,-6j-Dj 

with  a,  B,  N-i,...,  N^  specified.  Notice  that  the  optimal  6^  is  independent  of  a and  s. 
To  measure  the  extent  of  variance  reduction  with  a pe'^ticular  model, 

^U)  ^ V^'^^/n  min  [var^'^^(R^)] 


one  can  use 
(14) 


{ k) 

The  quantity  ni_^  ' denotes  the  number  of  replications  that  independent 
sampling  requires  to  achieve  the  same  variance  that  results  from  using  the 
optimal  sampling  plan  for  model  k.  For  n«de1  1 one  has 

(15)  u<’>  = 1/  |l  ♦ 2 (a-c)^  [C'"-  N,  " N,]  /V^’H 

; " i=1  J " j 

which  for  large  Np...,  N^  is 


(1 ) / , \ P n-1  'I 

(16)  u 1/[1  + 2(ci'^-  : N,  ;■; 


i=l 


N.)/p(l-p)  E NT] 
j=i+l  J i=l  ’ 


■ (l-p)/(l-np  + /N^np) 


where  the  upper  bound  obtains  when  = . . . = - N . 

Since  large  lead  to  a large  variance  for  when  using 


independent  replications,  we  explore  this  case  in  more  detail.  In  particular 
this  allows  us  to  ignore  linear  terms  in  in  ( ).  Supp^^se  that  p < 1/n. 


If  0^  (i'l)o  , where  p >_  o ^ 1/n  for  i = l 
so  that  ^ = 0 and 


n,  then 


cov  (I., 


♦ * * • > 
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(.y.  n-1  n n 5 

07)  -v.  l/[  l-2p  E N.  E N./O-P)  r Nf] 

" i=l  ’ j=i+l  J i=l  ^ 

= (1-P)/[1-P(  " N.)^/  E Nf]  < (l-p)/(l-np)  . 
i=l  ’ i=l  ^ 

This  case  avoids  the  need  to  solve  the  linear  program  in  (13  ) and  ( 14  ). 
Moreover  it  applies  equally  well  for  p i 1 - 1/n  since  one  can  substitute 
(1-p),  e and  o for  p,  o and  respectively, without  loss  of  generality. 

Solutions  of  the  linear  program  for  1/n  < p < 1/2  with  = ...  * N 

also  yield  worthwhile  insights.  Perusal  of  Table  1 for 

i/n  < p < (i+l)/n,  i = 1,‘...,  LP/2-lJ  showed  that  for  large  N 

(18)  var^^^R^)  = (a-e)^  q(l-nq)/n 

for  q = p-i/n  so  that 

(19)  = p(l-p)/q(l-nq)  . 

Tl-  variance  in  (18  ) vanishes  for  q=0  and  1/n  and  has  a maximum 
(a-3)^/4n^  for  q = l/2n  . This  implies  that 

(20)  min  var^^^  (R  ) < N^(a-3)^  /4n^  0 < p < 1/2 

i 

so  that 

(21 ) 0)^^ 4n  p(l-p)  . 

For  examples,  n=2  and  p=.25  leads  to  whereas  n=10  and  p=.4 

leads  to  ^^9.6  . 
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The  linear  program  solution  for  N^  = ...  = N^  = N also  reveals  that  the 
optimal  o_  need  not  have  distinct  entries.  For  example,  n = 8 and  p = .45 

gives  0-j  = 03  ” ®3  “ ®4  “ °5  ~ ®6  ” ®7  ~ ®8  ~ 

Here  one  forms  Rg  as 

(22)  Rq  = (2R^  + R3  + R4  + R5  2Rg  + Rg)/8. 


With  regard  to  the  actual  simulation,  one  would  use  Rg  = R^  and  R^  = Rg. 
These  results  imply  that  the  suggested  variance  reduction  technique  can  use 
no  more  than  six  replications  for  opt'*mality.  However,  the  weighting 
opportunities  when  n = 8 apparently  allow  a greater  variance  reduction  8.25 
than  in  the  case  of  n = 6,  7.07.  Some  additional  cases,  such  as  n = 9 and 
p = .20  also  need  clarification.  Notice  that  = 0 and  Og  = 1.  Since 
W9  . U + 0g  (mod  1)  we  have 

(23)  ITg  = (2R^  + 2Rg  + 2R^  + 2Rg  + Rg)/9 


Here  five  replications  suffice  for  optimality. 

Let  us  now  sumtiarize  the  properties  of  model  1 . 

1.  The  optimal  £ is  independent  of  a and  s.  See  (12). 

2.  For  large  N^,...,  is  independent  of  o and  B.  See  (16). 

3.  For  given  p increases  as  N^  ->■  N for  i = 1,...,  n.  See  (16) 

4.  For  equal  N^.  ~ ‘ '’'l)  where 

q r p (mod  1/n).  See  Table  1. 

★ 

5.  For  equal  the  optimal  solution  may  call  for  n < n repli- 
cations but  with  different  weights  for  each  observation.  The 
result  is  a smaller  variance  and  larger  variance  reduction  than 

it 

would  occur  for  the  optimal  solution  would  yield  for  n 
replications  with  equal  observation  weights.  See  Table  1. 
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3.  Model  2 

The  assumptions  in  model  1 that  X^...,  X^  be  i.i.d.,  Y^,...,  be  i.i.d. 
and  that  {X^. } and  {Y^ } be  independent  are  unnecessary  restrictions  whose  re- 
moval can  lead  to  an  additional  variance  reduction.  Assume  that  is 

« 

sufficiently  large  for  i=l,...,  n so  that  treating  X^  and  Y^.  as  normal 
variates  with  means  N^o  and  N^.6,  respectively,  and  variances  N^a(l-a)  and 
N^h(1-3),  respectively,  introduces  incidental  error.  Also  assume  that 
Wp...,  Wj.  follow  (5)  . Let  denote  a random  variable  from  the  normal 
distribution  with  zero  mean  and  unit  variance,  denoted  by  W(0,1).  Since 
either  X^.  or  Y^  occur  one  can  represent  these  quantities  by 


Xi  = N^a  + 


(24) 

Y.  = N.e  + z.  ^,.3(1-3)  . 
Then  (3)  recomes  for  replication  i 


aR^  = (a-3)N^  Al^  + A ^ Z^  Al^.  + (Ap+B)  /TiT  Z^. 


A = /o(l”a) 


B = /b(1-3^ 


- B 


so  that  has 


(25) 


var 


^^^(Rn)  = + 2ul^h/n 


where  as  is  expected,  and 

n n 


(26) 


(2). 


n-1  n 
y.  y 
i=l  j=i+l 


cov  (Ip  Ij)  [n.  H.  (a-3)^ 

+ A^  .liijliT  corr  (Zp  Z^)  ] * (Ap+B) ^.^5^  corr  (ZpZj)j 
If  Zp...,  Z^  are  independent  then 


(27) 
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(2)  2 
= I N.  z N cov  (I.,  I.)  (a-B)^ 
" i=l  ^j=i+l  ^ ^ 


However,  it  is  possible  to  n«ke  Z^,...,  Z^  negatively  correlated  so  ti’iat 
corr  (Zj,  Zj)  < 0 for  all  i^j.  In  particular,  the  most  negative  correlation 


achievable  is^ 


(28) 
so  that 
(29) 


ij^j 


(2) 

nU;'^^  = z z 
" i=l  3=i+l 


cov  (Z.,  1.)  = - l/(n-l) 

cov  (I.,  Ij)  [(a-B)^  H.  - aVn^  /(n-1)] 
- (Ap  + B)^ 

To  induce  (28  ) one  proceeds  as  follows.  Let  Zj,...,  1'^  be  independent 
from  W(0,  1).  Then  form 
(30) 

where  one  restricts  the  lower  triangular  matrix 


i=l ,. . . , n 


so  that 


(31, 


f h^^  0 ....  0 

!?12  t*22  ! 

• • • 

• * • , 

^ n ^2  n*‘*^nn 


z h^ij  = 1 


j=i 


Nk  'jk  = 


-l/(n-l) 


1 < J , i“l,...*  n 


In  particular,  solving  (31  ) yields 


(32) 


h^l  = 1 

h^j=-l/(n-l) 

j = 2 

" ^+l,j 

j “ 2 , . . . i j 

'''  This  is  seen  by  noting  that  any  smaller  covariance  makes  the  covariance 
matrix  of  Z^,...,  negative  definite . 
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Algorithtn  VRN  computes  h. 


Algorithm  VRfl(n) 

].  p <•  -l/(n-l)  . 

2.  1 . 

3.  For  i=2(l  )n  ■*-  p . 

4.  j ^ 2. 

5.  s 0. 

6.  For  i=l(l)j-l  s ^ s + h?. 

J I • 

7.  h..  --.'TT  . 

J J 

8.  If  j = n,  deliver  h. 

9.  s •♦-  p . 

10.  For  i = 1(l)j-l  s-H  s-h?. 

J * • 

12.  j -<•  j + 1 . 

13.  If  j < n,  for  i = j+l(l)n  h.  . , h.  . , 

' > J“  * J > J* I • 

14.  Go  to  5. 

Consider  the  case  for  which  cov(I^,  Ij)  = 0 i j.  Then 

(33)  = i/[i  - 2 (Ap  + B)^  (n  - 1)'^  E ^ 

” i=l  ^ j=i+l  J " 

is  the  variance  reduction.  Consider  the  case  N-j  = ...  = = N.  Then  (53)  is 


(34) 


= 1/n  - (Ap  + 0^/[p0  - P)(a  - e)^N  + pad  - a 


Expressions  (33)  and  (34)  allow  the  observations: 

1.  For  given  a ,Pand  p,  variance  reduction  decreases  as 

2.  For  p close  to  one  or  zero 

n n , n-1  n 


U) 


(2) 


i 


E N,  /[  E N.  - 2(n  - 1)”^  E E rfj. 

=1  ^ i=l  ^ i=l  ^ j=i+l  ^ 


) + (1  - p)3(l  - 3)]] 
Nj,...,  ‘1^  increase. 
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Clearly  this  model  is  of  most  benefit  when  the  mixture  is  heavily  .weighted 
toward  o or  p » a result  that  is  the  converse  of  an  observation  for  model  1 
where  small  p leads  to  sn*all  variance  reductions. 

4.  Hodel  3 

The  assumption  of  only  two  possible  survival  probabilities  a and  3 in  ? 
given  year  is  an  abstraction  usually  forced  on  an  investigator  by  incomplete 
knowledge.  In  a more  realistic  model  one  can  conceive  of  the  survival  probability 
P being  a continuous  function  of,  say,  rainfall  F.  For  example,  suppose  that 

(35)  P = ^ ^ a > 0,  b < 0 

where  F is  from  W(p,o)  and  pr(F  >^  - b/a)  0.  Here  ? increases  wHh  F 
as  one  might  expect  in  an  animal  population  where  rainfall  is  positively 
correlated  with  the  nutrUional  base  available  to  the  population. 

Let  subscripts  on  P and  F denote  replication  number  and  assume  that 

F^.  = 11  + Z^.o,  where  Z.  is  from  W(0,1).  If  Z is  formed  as  in  (30)  then 
corr  (F^,  F.)  = - l/(n-l;  for  n replications  and  one  can  show  that 

* V 

E(P^)  = f(l.l) 

f(4,2)  i = j 

f(2-2/{n-l),2)  i/j 

Let  denote  the  number  of  survivors  out  of  a population  of  size 
with  a randomly  selected  survival  probability  P^..  Then 


E(P,Pj)  - 
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(37) 


fE(R,|Pl)  = N^  ^ 

E(R^ |P^)  » N^  P^.(l  - p^)  + (N^  p^j2 

E(R|  RjlP^.Pj)  » N^  Nj.  P^  P^ 
E(R^)=f(l.l) 

E(R^)  = N.[f(M)  - f(4,2)  + N^f(4,2)] 
E(R^  Rj)  = N^  N.  f(2  - 2/(n  - 1).  2) 


i f j 


i f j 


so  that 


n 


(38)  n^  var^^^(Rj  = [f(4,2)  - f(2,2)]  £ N?  + [f(l ,1 ) - f(4,2)]  z H. 

" ;=]  ^ i=l  ’ 

+ 2[f(2  - 2/(n  - 1).  2)  - f(2,2)]  Z Z H.  H.. 

i=l  j=i+l  ’ J 


Since 


a 

* 


fc 


(39)  f(2  - 2/(n  - 1).  2)  = f(2,  2)f(-  2/(n  - 1).  0) 

the  covariance  contribution  in  (38)  is  negative. 

Consider  the  case  Ni  = ...  N * N for  which  the  variance  reduction  is 

I n 

(40)  I/O  + (n-l)[f(-2/(n-l),0)  - l]/[f(2,0)  -1  + (f(^l,-l)  -f(2.0))/N]}.- 
for  N large  this  is 

(41)  1/(1  + (n-l)[f(-2/(n-l),  0)  - l]/[f(2,0)  - 1]}  . 

2 

For  given  a and  a (41)  decreases  monotonically  as  n increases,  suggesting 
the  appeal  of  correlating  replications  in  pairs.  Moreover,  variance  considera- 
tions indicate  that  for  a total  of  n replications  creating  n/2  (n  even) 
independent  sets  of  replications  each  with  two  replications  correlated  according 
to  model  3 yields  a smaller  variance  than  using  model  3 to  induce  correlation 

among  all  n replications.  Also,  for  a given  n (41)  monotonically  decreases 

2 (31 

as  a and  o increase.  Table  2 shows  selected  values  of  • Notice  that 

the  greatest  variance  reduction  accrues  when  P^  varies  little  between  successive 

years. 


S.  Example^ 

By  way  of  illustration  we  describe  the  application  of  model  1 to  a Simula- 
tion  of  the  population  dynamics  of  a hypothetical  female  elephant  population. 
Three  attributes  characterize  each  strata:  age,  maturity  and  pregnancy.  Only 
mature  elephant  conceive  and  the  gestation  period  is  twenty- two  months.  Figure 
1 shows  how  strata  change  from  year  to  year.  Each  arc  has  an  associated  prob- 
ability. In  some  cases  this  probability  may  be  zero.  For  example,  one  year  old 
calves  have  zero  probability  of  maturing.  Other  probabilities  are  unity,  as  in 
the  case  of  surviving  pregnant  elephants  giving  birth. 

The  remaining  probabilities  are  functions  of  rainfall  and  density.  In 

particular,  survival,  maturity  and  conception  probabilities  increase  as  rainfall 
increases  but  decrease  as  population  density  increases.  For  expository  purposes 

the  example  neglects  density  dependence.  However  it  does  specify  explicit  re- 
lationships between  the  probabilities  and  rainfall.  The  variance  reduction 
technique  of  model  1 was  applied  to  the  survival  probabilities.  Specifically 

^ This  model  emerged  from  discussions  with  Professor  Daniel  Botkin  at  the 
Ecosystems  Center  at  the  Marine  Biological  Laboratory  at  Woods  Hole  and  wit;*- 
Professors  Richard  Miller  and  Matthew  Sobel  of  Yale  University. 


-18- 


let  F denote  rainfall  in  a given  year  and  assume  it  has  a normal  distribution. 

Then  the  survival  probability  is  Oj^  for  strata  k if  F ^ F*  and  if 

* * 

F < F , where  oj^  > 3|^.  The  quantity  F is  a threshold  rainfall  that  deter- 
mines survival  prospects.  In  terms  of  model  1 p = pr(F  > F*).  In  the 
present  example  p = 0.25. 

At  each  point  in  time  we  performed  10  independent  sets  of  sampling  experi- 
ments for  each  strata*  where  each  experiment  consisted  of  n=4  replications. 

The  independent  sets  enabled  us  to  estimate  the  variance  of  the  average  over 
the  4 replications  subjecte*'  to  the  variance  reduction  technique  in  each  strata. 
For  p = .25,  0.  = 0,  02  = .25,  0^  = .50  and  0^  = .75  are  an  optimal  solution. 
Each  simulation  began  with  the  same  population  profile  and  ran  for  100  years. 

For  a given  strata  define  R..  as  the  number  of  surviving  elephants  on  replica- 

“ *} 

tion  i of  set  j.  Then 


(43) 
where 

(44) 


1 

40  X 3$ 


4 

z 

i=l 


10 

l (R 
j=l 


IJ 


- R..)^ 


R 


1 

10 


10 

z 

j=l 


» 


provides  an  estimate  of  the  variance  of  R..  if  all  replications  were  independent 
and 

(45)  T = ygjy  Z (R.j  - R..)^ 

V * 

provides  an  estimate  of  the  variance  of  R..  regardless  of  whether  or  not  the 
replications  are  correlated.  Then  for  strata  k Sj^/Tj^  provides  a measure  of 
variance  reduction.  Since  sampling  in  strata  k independent  of  sampling 
in  strata  £ (k  / £)  in  a given  year  a summary  measure  of  variance  reduction  is 
(>:S|^)/(zTj^).^  Table  2 shows  this  ratio  for  years  m = 5(5)100.  In  general, 


1 The  only  contrary  case  occurs  between  pregnant  elephants  that  give  birth  and 
the  new  born  calves.  However,  neglecting  this  correlation  should  not  be  a 
serious  issue. 
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